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Abstract 

In this paper, we review the construction of low-dimensional manifolds of reduced 
description for equations of chemical kinetics from the standpoint of the method 
of invariant manifold (MIM). MIM is based on a formulation of the condition of 
invariance as an equation, and its solution by Newton iterations. A review of exist- 
ing alternative methods is extended by a thermodynamically consistent version of 
the method of intrinsic low-dimensional manifolds. A grid-based version of MIM is 
developed, and model extensions of low-dimensional dynamics are described. Gen- 
eralizations to open systems are suggested. The set of methods covered makes it 
possible to effectively reduce description in chemical kinetics. 

Key words: Kinetics, Model Reduction, Invariant Manifold, Entropy, Nonlinear 
Dynamics, Mathematical Modeling 



1 Introduction 



In this paper, we present a general method of constructing the reduced de- 
scription for dissipative systems of reaction kinetics. Our approach is based on 
the method of invariant manifold which was developed in end of 1980th - be- 
ginning of 1990th by |Gorban fc Karlin (199^ fl[b|). Its realization for a generic 
dissipative systems was discussed by |Gorban fc Karlin (1994] ) ; |Gorban, Kar 



[lin, Ilg fc Ottinger (2001] ). This method was applied to a set of problems of 
classical kinetic theory based on the Boltzmann kinetic equation (porbarTS; 
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Karlin, 1994; [Karlin, Dukek fc Nonnenmacher, 1997] ; [Karlin, Gorban, Dukek| 



fc Nonncnmacher, 1998 ). The method of invariant manifold was successfully 
applied to a derivation of reduced description for kinetic equations of poly- 
meric solutions flZmievskii, Kalin fc Deville, 2000|) . It was also been tested on 



systems of chemical kinetics ( |Gorban, Karlin, Zmievskii fc Dymova, 2000| ). 

The goal of nonequilibrium statistical physics is the understanding of how a 
system with many degrees of freedom acquires a description with a few degrees 
of freedom. This should lead to reliable methods of extracting the macroscopic 
description from a detailed microscopic description. 

Meanwhile this general problem is still far from the final solution, it is reason- 
able to study simplified models, where, on the one hand, a detailed description 
is accessible to numerics, on the other hand, analytical methods designed to 
the solution of problems in real systems can be tested. 

In this paper we address the well known class of finite-dimensional systems 
known from the theory of reaction kinetics. These are equations governing 
a complex relaxation in perfectly stirred closed chemically active mixtures. 
Dissipative properties of such systems are characterized with a global con- 
vex Lyapunov function G (thermodynamic potential) which implements the 
second law of thermodynamics: As the time t tends to infinity, the system 
reaches the unique equilibrium state while in the course of the transition the 
Lyapunov function decreases monotonically. 

While the limiting behavior of the dissipative systems just described is cer- 
tainly very simple, there are still interesting questions to be asked about. One 
of these questions is closely related to the above general problem of nonequi- 
librium statistical physics. Indeed, evidence of numerical integration of such 
systems often demonstrates that the relaxation has a certain geometrical struc- 
ture in the phase space. Namely, typical individual trajectories tend to man- 
ifolds of lower dimension, and further proceed to the equilibrium essentially 
along these manifolds. Thus, such systems demonstrate a dimensional reduc- 
tion, and therefore establish a more macroscopic description after some time 
since the beginning of the relaxation. 



There are two intuitive ideas behind our approach, and we shall now discuss 
them informally. Objects to be considered below are manifolds (surfaces) f2 
in the phase space of the reaction kinetic system (the phase space is usually a 
convex polytope in a finite-dimensional real space). The 'ideal' picture of the 
reduced description we have in mind is as follows: A typical phase trajectory, 
c(t), where t is the time, and c is an element of the phase space, consists of 
two pronounced segments. The first segment connects the beginning of the 
trajectory, c(0), with a certain point, c(ti), on the manifold Q (rigorously 
speaking, we should think of c(ti) not on f2 but in a small neighborhood of 
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n but this is inessential for the ideal picture). The second segment belongs to 
f2, and connects the point c(t\) with the equilibrium c eq = c(oo), c eq e f2. 
Thus, the manifolds appearing in our ideal picture are "patterns" formed by 
the segments of individual trajectories, and the goal of the reduced description 
is to "filter out" this manifold. 

There are two important features behind this ideal picture. The first feature is 
the invariance of the manifold f2: Once the individual trajectory has started 
on ft, it does not leaves f2 anymore. The second feature is the projecting: The 
phase points outside fl will be projected onto Q. Furthermore, the dissipativity 
of the system provides an additional information about this ideal picture: 
Regardless of what happens on the manifold f2, the function G was decreasing 
along each individual trajectory before it reached f2. This ideal picture is the 
guide to extract slow invariant manifolds. 

The paper is organized as follows. In the section 2, we review the reaction 
kinetics (section 2.1), and discuss the main methods of model reduction in 
chemical kinetics (section 2.2). In particular, we present two general versions 
of extending partially equilibrium manifolds to a single relaxation time model 
in the whole phase space, and develop a thermodynamically consistent version 
of the intrinsic low-dimensional manifold (ILDM) approach. In the section 3 
we introduce the method of invariant manifold in the way appropriate to this 
class of nonequilibrium systems. In the sections 4 and 5 we give some details 
on the two relatively independent parts of the method, the thermodynamic 
projector, and the iterations for solving the invariance equation. We also intro- 
duce a general symmetric linearization procedure for the invariance equation, 
and discuss its relevance to the picture of decomposition of motions. In the 
section 6, these two procedures are combined into an unique algorithm. In the 
section 7, we demonstrate an illustrative example of computations for a model 
catalytic reaction. In the section 8 we demonstrate how the thermodynamic 
projector is constructed without the a priori parameterization of the manifold. 
This result is essentially used in the section 9 where we introduce a computa- 
tionally effective grid-based method to construct invariant manifolds. In the 
section 10 we describe an extension of the method of invariant manifold to 
open systems. Finally, results are discussed in the section 11. 



2 Equations of chemical kinetics and their reduction 

2.1 Outline of the dissipative reaction kinetics 



We begin with an outline of the reaction kinetics (for details see e. g. the book 
of |Yablonskii, Bykov, Gorban fc Elokhin (19"9~T| )). Let us consider a closed sys- 
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tern with n chemical species Ai, . . . , A n , participating in a complex reaction. 
The complex reaction is represented by the following stoichiometric mecha- 
nism: 

a s iAi + . . . + a sn A n ^ /3 sl Ai + ... + (3 sn k n , (1) 

where the index s = 1, . . . ,r enumerates the reaction steps, and where in- 
tegers, a S i and f3 S i, are stoichiometric coefficients. For each reaction step s, 
we introduce n-component vectors ct s and (3 S with components a S i and (3 S i. 
Notation 7 S stands for the vector with integer components 7 S j = f3 S i — a S i (the 
stoichiometric vector). We adopt an abbreviated notation for the standard 
scalar product of the n-component vectors: 

n 

The system is described by the n-component concentration vector c, where the 
component q > represents the concentration of the specie Aj. Conservation 
laws impose linear constraints on admissible vectors c (balances): 

(bi,c) = Bi, i = I,..., I, (2) 

where bi are fixed and linearly independent vectors, and Bi are given scalars. 
Let us denote as B the set of vectors which satisfy the conservation laws (2): 

B = {c\(b 1 ,c} = B 1 ,...,(b l ,c)=B l }. 

The phase space V of the system is the intersection of the cone of n-dimensional 
vectors with nonnegative components, with the set B, and dimV — d — n — l. 
In the sequel, we term a vector c e V the state of the system. In addition, 
we assume that each of the conservation laws is supported by each elementary 
reaction step, that is 

(7 S A>=0, (3) 

for each pair of vectors 7 S and 6j. 

Reaction kinetic equations describe variations of the states in time. Given the 
stoichiometric mechanism (1), the reaction kinetic equations read: 

c = J(c), J(c) = 5>.W.(c), (4) 

s=l 



4 



where dot denotes the time derivative, and W s is the reaction rate function of 
the step s. In particular, the mass action law suggests the polynomial form of 
the reaction rates: 

n n 

w. = Kn<?-k:i[<*, (5) 

1=1 1=1 

where kf and k~ are the constants of the direct and of the inverse reactions 
rates of the sth reaction step. The phase space V is positive-invariant of the 
system (4): If c(0) G V, then c(t) G V for all the times t > 0. 

In the sequel, we assume that the kinetic equation (4) describes evolution 
towards the unique equilibrium state, c eq , in the interior of the phase space 
V. Furthermore, we assume that there exists a strictly convex function G(c) 
which decreases monotonically in time due to Eq. (4): 

G=(VG(c),J(c))<0, (6) 

Here VG is the vector of partial derivatives dGj dQ, and the convexity assumes 
that the n x n matrices 

H c =\\&G(c)/dc i dc j \\, (7) 

are positive definite for all c G V. In addition, we assume that the matrices 
(7) are invertible if c is taken in the interior of the phase space. 

The function G is the Lyapunov function of the system (4), and c cq is the 
point of global minimum of the function G in the phase space V. Otherwise 
stated, the manifold of equilibrium states c eq (f?i, . . . , B{) is the solution to the 
variational problem, 

G — > min for (b i: c) — B i: i = 1, . . . , /. (8) 

For each fixed value of the conserved quantities B h the solution is unique. 
In many cases, however, it is convenient to consider the whole equilibrium 
manifold, keeping the conserved quantities as parameters. 

For example, for perfect systems in a constant volume under a constant tem- 
perature, the Lyapunov function G reads: 

n 

G = E^NQ/cD-l]. (9) 
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It is important to stress that c eq in Eq. (9) is an arbitrary equilibrium of the 
system, under arbitrary values of the balances. In order to compute G(c), 
it is unnecessary to calculate the specific equilibrium c cq which corresponds 
to the initial state c. Moreover, for ideal systems, function G is constructed 
from the thermodynamic data of individual species, and, as the result of this 
construction, it turns out that it has the form of Eq. (9). Let us mention here 
the classical formula for the free energy F = RTVG: 

F = VRTj2ci[(HciV Q - 1) + F int i(T)], (10) 

i=l 



where V is the volume of the system, T is the temperature, Vq j = No(2wh 2 /mikT) 3 / 2 
is the quantum volume of one mole of the specie Aj, Nq is the Avogadro num- 
ber, rrii is the mass of the molecule of Aj, R = kN , and -F int j(T) is the free 
energy of the internal degrees of freedom per mole of Aj. 

Finally, we recall an important generalization of the mass action law (5), 
known as the Marcelin-De Donder kinetic function. This generalization was 
developed by [Feinberg [1972) based on ideas of the thermodynamic theory of 



affinity ( pe Donder fc Van Kysselberghe, 193(j| ). We use the kinetic function 



suggested by |Bykov, Gorban fc Yablonskii (19821) . Within this approach, the 



functions W s are constructed as follows: For a given strictly convex function 
G, and for a given stoichiometric mechanism (1), we define the gain (+) and 
the loss (— ) rates of the sth step, 

^ + = y?+exp[(VG> s >], W" = <pj exp[(VG, /3 S >], (11) 



where (pf > are kinetic factors. The Marcelin-De Donder kinetic function 
reads: W s = — W~, and the right hand side of the kinetic equation (4) 
becomes, 

J = jl IsWt exp[( VG, a,)] - <p~ exp[(VG, &>]}. (12) 

8=1 



For the Marcelin-De Donder reaction rate (11), the dissipation inequality (6) 
reads: 

G = ±[(VG,(3 S ) - (VG,a s >] Lfe^^-^e^A < 0. (13) 

s=l ^ J 



The kinetic factors ipf should satisfy certain conditions in order to make valid 
the dissipation inequality (13). A well known sufficient condition is the detail 
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balance: 



(14) 



other sufficient conditions are discussed in detail elsewhere flYablonskii, Bykov 
Gorban fc Elokhin, 19911 ; |Gorban, 198j |Karlin, 1989j , p93|) . For the function 
G of the form (9), the Marcelin-De Donder equation casts into the more fa- 
miliar mass action law form (5). 



2.2 The problem of reduced description in chemical kinetics 



What does it mean, "to reduce the description of a chemical system"? This 



means the following: 



(1) 



(2) 



(3) 



To shorten the list of species. This, in turn, can be achieved in two ways: 

(i) To eliminate inessential components from the list; 

(ii) To lump some of the species into integrated components. 

To shorten the list of reactions. This also can be done in several ways: 

(i) To eliminate inessential reactions, those which do not significantly 
influence the reaction process; 

(ii) To assume that some of the reactions "have been already com- 
pleted", and that the equilibrium has been reached along their paths 
(this leads to dimensional reduction because the rate constants of the 
"completed" reactions are not used thereafter, what one needs are equi- 
librium constants only). 

To decompose the motions into fast and slow, into independent (almost- 
independent) and slaved etc. As the result of such a decomposition, the 
system admits a study "in parts". After that, results of this study are 
combined into a joint picture. There are several approaches which fall 
into this category: The famous method of the quasi-steady state (QSS), 
pioneered by Bodenstein and Semenov and explored in considerable detail 
by many authors, in particular, by |Bowen, Acrivos fc Oppenheim (1963|) ; 
|Ghen (1988| ); frjegel k Slemrod (1989| ); |Fraser (1988| ); |Roussel fc Fragei 
(1990| ), and many others; the quasi-equilibrium approximation ( Prlov & 



Rozonoer, 1981 IGorban, 1981 IVolpert fc Hud.jaev, 19851 ; IFraser, 1988) ; 



Karlin, 1989 



methods of sensitivity analysis ( [Rabitz, Kramer"^ 



Dacol, 1983| ; [Lam fc Goussis, 1994Q ; methods based on the derivation of 



the so-called intrinsic low- dimensional manifolds (ILDM, as suggested by 
[Vlaas fc Pope (1992| )). Our method of invariant manifold (MIM, ( |Gorban| 
fc Karlin, 1992| Ja|Jb], |1994| ; |Gorban, Karlin, Zmievskii fc Dymova, 200C| ; 



porban, Karlin, Ilg fc Ottinger, 2001| )) also belongs to this kind of meth- 
ods. 
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Why to reduce description in the times of supercomputers? 



First, in order to gain understanding. In the process of reducing the description 
one is often able to extract the essential, and the mechanisms of the processes 
under study become more transparent. Second, if one is given the detailed 
description of the system, then one should be able also to solve the initial- 
value problem for this system. But what should one do in the case where the 
the system is representing just a point in a three-dimensional flow? The prob- 
lem of reduction becomes particularly important for modeling the spatially 
distributed physical and chemical processes. Third, without reducing the ki- 
netic model, it is impossible to construct this model. This statement seems 
paradoxal only at the first glance: How come, the model is first simplified, 
and is constructed only after the simplification is done? However, in practice, 
the typical for a mathematician statement of the problem, (Let the system 
of differential equations be given, then ...) is rather rarely applicable in the 
chemical engineering science for detailed kinetics. Some reactions are known 
precisely, some other - only hypothetically. Some intermediate species are well 
studied, some others - not, it is not known much about them. Situation is 
even worse with the reaction rates. Quite on the contrary, the thermodynamic 
data (energies, enthalpies, entropies, chemical potentials etc) for sufficiently 
rarefied systems are quite reliable. Final identification of the model is always 
done on the basis of comparison with the experiment and with a help of fit- 
ting. For this purpose, it is extremely important to reduce the dimension of the 
system, and to reduce the number of tunable parameters. The normal logics 
of modeling for the purpose of chemical engineering science is the following: 
Exceedingly detailed but coarse with respect to parameters system — > reduc- 
tion — > fitting — > reduced model with specified parameters (cycles are allowed 
in this scheme, with returns from fitting to more detailed models etc). A more 
radical viewpoint is also possible: In the chemical engineering science, detailed 
kinetics is impossible, useless, and it does not exist. For a recently published 
discussion on this topic see [Levenspiel (1999| , |2000| ) ; [Yablonsky (2000|) . 



Alas, with a mathematical statement of the problem related to reduction, we 
all have to begin with the usual: Let the system of differential equations be 
given .... Enormous difficulties related to the question of how well the original 
system is modeling the real kinetics remain out of focus of these studies. 



Our present work is devoted to studying reductions in a given system of ki- 
netic equations to invariant manifolds of slow motions. We begin with a brief 
discussion of existing approaches. 
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2.3 Partial equilibrium approximations 

Quasi- equilibrium with respect to reactions is constructed as follows: From the 
list of reactions (1), one selects those which are assumed to equilibrate first. Let 
they be indexed with the numbers s±, . . . , s^. The quasi-equilibrium manifold 
is defined by the system of equations, 

W+ = W~, i = l,...,k. (15) 

This system of equations looks particularly elegant when written in terms of 
conjugated (dual) variables, /x = VG: 

( 7si ,/x)=0, i — l,...,k. (16) 

In terms of conjugated variables, the quasi-equilibrium manifold forms a linear 
subspace. This subspace, L- 1 , is the orthogonal completement to the linear 
envelope of vectors, L = lin{7 si , . . . , 7 Sfc }- 

Quasi-equilibrium with respect to species is constructed practically in the same 
way but without selecting the subset of reactions. For a given set of species, 
A il , . . . , A ik , one assumes that they evolve fast to equilibrium, and remain 
there. Formally, this means that in the /c-dimensional subspace of the space of 
concentrations with the coordinates , . . . , q fc , one constructs the subspace L 
which is defined by the balance equations, (&«, c) = 0. In terms of the conju- 
gated variables, the quasi-equilibrium manifold, L- 1 , is defined by equations, 

H e L x , (/i= (//!,...,//„)). (17) 

The same quasi-equilibrium manifold can be also defined with the help of 
fictitious reactions: Let g 1 ,...,g be a basis in L. Then Eq. (17) may be 
rewritten as follows: 

(^,Ai> = 0, i = l,...,q. (18) 

Illustration: Quasi-equilibrium with respect to reactions in hydrogen oxidation: 
Let us assume equilibrium with respect to dissociation reactions, H 2 ^ 2H, 
and, O2 ^ 20, in some subdomain of reaction conditions. This gives: 

kt c H 2 = ki c H , /cJco 2 = k 2 Cq. 

Quasi-equilibrium with respect to species: For the same reaction, let us assume 
equilibrium over H, O, OH, and H 2 2 , in a subgomain of reaction conditions. 
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Subspace L is defined by balance constraints: 



ch + coh + 2c H2 o 2 = 0, c + c H + 2c H2 o 2 = 0. 

Subspace L is two-dimensional. Its basis, {gi,g 2 } in the coordinates ch, Co, 
coh, and ch 2 o 2 reads: 

9l = (1,1,-1,0), g 2 = (2,2,0,-1). 

Corresponding Eq. (18) is: 

A*H + A*0 = /^OH, 2/i.H + 2/i.Q = A i H 2 o 2 - 



General construction of the quasi- equilibrium manifold: In the space of con- 
centration, one defines a subspace L which satisfies the balance constraints: 

(bi,L) = 0. 



The orthogonal complement of L in the space with coordinates n = VG 
defines then the quasi-equilibrium manifold Ql. For the actual computations, 
one requires the inversion from fi to c. Duality structure ii <-> c is well studied 
by many authors ( prlov fc Rozonoer, 1984] ; pukek, Karlin fc JNonnenmacher 

"15571 ). 



Quasi-equilibrium projector. It is not sufficient to just derive the manifold, it 
is also required to define a projector which would transform the vector field 
defined on the space of concentrations to a vector field on the manifold. Quasi- 
equilibrium manifold consists of points which minimize G on the affine spaces 
of the form c+L. These affine planes are hypothetic planes of fast motions (G is 
decreasing in the course of the fast motions). Therefore, the quasi-equilibrium 
projector maps the whole space of concentrations on parallel to L. The 
vector field is also projected onto the tangent space of Ql parallel to L. 

Thus, the quasi-equilibrium approximation implies the decomposition of mo- 
tions into the fast - parallel to L, and the slow - along the quasi-equilibrium 
manifold. In order to construct the quasi-equilibrium approximation, knowl- 
edge of reaction rate constants of "fast" reactions is not required (stoichio- 
metric vectors of all these fast reaction are in L, 7 fast G L, thus, knowledge 
of L suffices), one only needs some confidence in that they all are sufficiently 
fast ([Volpert fc Hudjaev, 1985|) . The quasi-equilibrium manifold itself is con- 



structed based on the knowledge of L and of G. Dynamics on the quasi- 
equilibrium manifold is defined as the quasi-equilibrium projection of the "slow 
component" of kinetic equations (4). 
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2.4 Model equations 



The assumption behind the quasi-equilibrium is the hypothesis of the decom- 
position of motions into fast and slow. The quasi-equilibrium approximation 
itself describes slow motions. However, sometimes it becomes necessary to re- 
store to the whole system, and to take into account the fast motions as well. 
With this, it is desirable to keep intact one of the important advantages of 
the quasi-equilibrium approximation - its independence of the rate constants 
of fast reactions. For this purpose, the detailed fast kinetics is replaced by a 
model equation (single relaxation time approximation) . 

Quasi-equilibrium models (QEM) are constructed as follows: For each concen- 
tration vector c, consider the affine manifold, c + L. Its intersection with the 
quasi-equilibrium manifold Ql consists of one point. This point delivers the 
minimum to G on c + L. Let us denote this point as c* L (c). The equation of 
the quasi-equilibrium model reads: 

c = ~\ [c- cl(c)] + £ l s W s {c* L {c)\ (19) 

slow 



where r > is the relaxation time of the fast subsystem. Rates of slow re- 
actions are computed in the points c* L (c) (the second term in the right hand 
side of Eq. (19), whereas the rapid motion is taken into account by a sim- 
ple relaxational term (the first term in the right hand side of Eq. (19). The 
most famous model kinetic equation is the BGK equation in the theory of 
the Boltzmann equation ( Bhatnagar, Gross fc Krook, 1954| ). The general the- 



ory of the quasi-equilibrium models, including proofs of their thermodynamic 
consistency, was constructed by |Gorban fc Karlin (1992c| , |1994a| ). 



Single relaxation time gradient models (SRTGM) were considered by [Ansumali 



|fc Karlin (2000| , |2002| ,|a|) in the context of the lattice Boltzmann method for 
hydrodynamics. These models are aimed at improving the obvious drawback 
of quasi-equilibrium models (19): In order to construct the QEM, one needs 
to compute the function, 

ct (c) = are min G(x). (20) 



This is a convex programming problem. It does not always has a closed-form 
solution. 

Let g 1 ,...,g k is the orthonormal basis of L. We denote as D(c) the k x 
k matrix with the elements {g { , Hc9j) , where He is the matrix of second 
derivatives of G (7). Let C(c) be the inverse of D(c). The single relaxation 
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time gradient model has the form: 

c= -\Y,9 i C{c) ij {g j , VG> + E 7s %). (21) 

i,j slow 

The first term drives the system to the minimum of G on c + L, it does not 
require solving the problem (20), and its spectrum in the quasi-equilibrium is 
the same as in the quasi-equilibrium model (19). Note that the slow component 
is evaluated in the "current" state c. 

The models (19) and (21) lift the quasi-equilibrium approximation to a kinetic 
equation by approximating the fast dynamics with a single "reaction rate 
constant" - relaxation time r. 



2.5 Quasi-steady state approximation 



The quasi-steady state approximation (QSS) is a tool used in a huge amount 
of works. Let us split the list of species in two groups: The basic and the 
intermediate (radicals etc). Concentration vectors are denoted accordingly, c s 
(slow, basic species), and c f (fast, intermediate species). The concentration 
vector c is the direct sum, c = c s © c f . The fast subsystem is Eq. (4) for the 
component c f at fixed values of c s . If it happens that this way defined fast 
subsystem relaxes to a stationary state, c f — > c qss (c s ), then the assumption 
that c f = Cq SS (c) is precisely the QSS assumption. The slow subsystem is the 
part of the system (4) for c s , in the right hand side of which the component 
c f is replaced with Cq SS (c). Thus, J = J s © Jf, where 



c f = J f (c s © c f ), c s = const; c f -> c { qss (c s ); (22) 
t s = J s (c*®c qss (c s )). (23) 

Bifurcations in the system (22) under variation of c s as a parameter are con- 
fronted to kinetic critical phenomena. Studies of more complicated dynamic 
phenomena in the fast subsystem (22) require various techniques of averaging, 
stability analysis of the averaged quantities etc. 

Various versions of the QSS method are well possible, and are actually used 
widely, for example, the hierarchical QSS method. There, one defines not a 
single fast subsystem but a hierarchy of them, c fl , . . . , c ffc . Each subsystem c fi 
is regarded as a slow system for all the foregoing subsystems, and it is regarded 
as a fast subsystem for the following members of the hierarchy. Instead of one 
system of equations (22), a hierarchy of systems of lower- dimensional equations 
is considered, each of these subsystem is easier to study analytically. 
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Theory of singularly perturbed systems of ordinary differential equations is 
used to provide a mathematical background and further development of the 
QSS approximation ( [Bowen, Acrivos fc Oppenheim, 1963] ; |Segel fc Slemrod 



1989| ). In spite of a broad literature on this subject, it remains, in general, un- 
clear, what is the smallness parameter that separates the intermediate (fast) 
species from the basic (slow). Reaction rate constants cannot be such a param- 
eter (unlike in the case of the quasi-equilibrium) . Indeed, intermediate species 
participate in the same reactions, as the basic species (for example, H 2 ^ 2H, 
H + O2 ^ OH + O). It is therefore incorrect to state that c f evolve faster than 
c s . In the sense of reaction rate constants, c f is not faster. 

For catalytic reactions, it is not difficult to figure out what is the smallness 
parameter that separates the intermediate species from the basic, and which 
allows to upgrade the QSS assumption to a singular perturbation theory rigor- 
ously flYablonskii. Bykov, Gorban fc Eiokhin, 1991 ). This smallness parameter 



is the ratio of balances: Intermediate species include the catalyst, and their 
total amount is simply significantly less than the amount of all the q's. After 
renormalizing to the variables of one order of magnitude, the small parameter 
appears explicitly. 

For usual radicals, the origin of the smallness parameter is quite similar. There 
are much less radicals than the basic species (otherwise, the QSS assumption is 
inapplicable). In the case of radicals, however, the smallness parameter cannot 
be extracted directly from balances Bi (2). Instead, one can come up with a 
thermodynamic estimate: Function G decreases in the course of reactions, 
whereupon we obtain the limiting estimate of concentrations of any specie: 

Cj < max c,-, (24) 

G(C)<G(C(0)) 



where c(0) is the initial composition. If the concentration cr of the radical 
R is small both initially and in the equilibrium, then it should remain small 
also along the path to the equilibrium. For example, in the case of ideal G (9) 
under relevant conditions, for any t > 0, the following inequality is valid: 

c R [ln(c R (t)/c e R q ) - 1] < G(c(0)). (25) 



Inequality (25) provides the simplest (but rather coarse) thermodynamic es- 
timate of cr(£) in terms of G(c(0)) and uniformly for t > 0. Complete 
theory of thermodynamic estimates of dynamics has been developed by |Gor- 



ban (1984] ). One can also do computations without a priori estimations, if one 
accepts the QSS assumption until the values c f stay sufficiently small. 

Let us assume that an a priori estimate has been found, Cj(t) < q max , for 
each Q. These estimate may depend on the initial conditions, thermodynamic 
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data etc. With these estimates, we are able to renormalize the variables in 
the kinetic equations (4) in such a way that renormalized variables take their 
values from the unit segment [0, 1]: q = q/cj max . Then the system (4) can be 
written as follows: 

§ = — 4(c). (26) 

CLZ Ci rnax 



The system of dimensionless parameters, e« = Cj max / maxj q max defines a 
hierarchy of relaxation times, and with its help one can establish various re- 
alizations of the QSS approximation. The simplest version is the standard 
QSS assumption: Parameters £j are separated in two groups, the smaller ones, 
and of the order 1. Accordingly, the concentration vector is split into c s © c f . 
Various hierarchical QSS are possible, with this, the problem becomes more 
tractable analytically. 

Corrections to the QSS approximation can be addressed in various ways (see, 
e. g., [Vasireva, Butuzov fc Kalachev (1995p ; |Strygin fc Sobolev (1988|) ). There 



exist a variety of ways to introduce the smallness parameter into kinetic equa- 
tions, and one can find applications to each of the realizations. However, the 
two particular realizations remain basic for chemical kinetics: (i) Fast reac- 
tions (under a given thermodynamic data); (ii) Small concentrations. In the 
first case, one is led to the quasi-equilibrium approximation, in the second 
case - to the classical QSS assumption. Both of these approximations allow 
for hierarchical realizations, those which include not just two but many relax- 
ation time scales. Such a multi-scale approach essentially simplifies analytical 
studies of the problem. 

The method of invariant manifold which we present below in the section 6 
allows to use both the QE and the QSS as initial approximations in the it- 
erational process of seeking slow invariant manifolds. It is also possible to 
use a different initial ansatz chosen by a physical intuition, like, for example, 
the Tamm-Mott-Smith approximation in the theory of strong shock waves 
QGorban fc Karlin, 1992Q . 



2.6 Methods based on spectral decomposition of Jacobian fields 



The idea to use the spectral decomposition of Jacobian fields in the problem of 
separating the motions into fast and slow originates from methods of analysis 
of stiff systems ( pear, 1971|) , and from methods of sensitivity analysis in con- 



trol theory ( [Rabitz, Kramer fc Dacol, 1983| ). There are two basic statements of 
the problem for these methods: (i) The problem of the slow manifold, and (ii) 
The problem of a complete decomposition (complete integrability) of kinetic 
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equations. The first of these problems consists in constructing the slow mani- 
fold f2, and a decomposition of motions into the fast one - towards f2, and the 
slow one - along f2 ( [Maas fc Pope, 1992| ). The second of these problems con- 
sists in a transformation of kinetic equations (4) to a diagonal form, Q = fi(Q) 
(so-called full nonlinear lumping or modes decoupling, [Lam fc Goussis (1994| ); 



[Li, Rabitz fc Toth (19941) ; Toth, Li, Rabitz fc Tomlin (1997|) ). Clearly, if one 



finds a sufficiently explicit solution to the second problem, then the system 
(4) is completely integrable, and nothing more is needed, the result has to 
be simply used. The question is only to what extend such a solution can be 
possible, and how difficult it would be as compared to the first problem to 
find it. 

One of the currently most popular methods is the construction of the so- 
called intrinsic low- dimensional manifold (ILDM, [Maas fc Pope (1992| )). This 



method is based on the following geometric picture: For each point c, one 
defines the Jacobian matrix of Eq. (4), Fc = dJ(c)/dc. One assumes that, in 
the domain of interest, the eigenvalues of Fc are separated into two groups, 
A^ and A!-, and that the following inequalities are valid: 



Re A • > a > 6 > ReAj, a > b, b < 0. 



Let us denote as L c and L c the invariant subspaces corresponding to X s and 
A f , respectively, and let Z s c and Z c be the corresponding spectral projectors, 
Z S C L C = L c , Z l c U c = L c , Z s c L f c = Z C L C = {0}, Z c + Z l c = 1. Operator Z c 
projects onto the subspace of "slow modes" L s c , and it annihilates the "fast 
modes" L c . Operator Z c does the opposite, it projects onto fast modes, and 
it annihilates the slow modes. The basic equation of the ILDM reads: 

Z f c J(c) = 0. (27) 



In this equation, the unknown is the concentration vector c. The set of solu- 
tions to Eq. (27) is the ILDM manifold fiiidm- 

For linear systems, F c , Z s c , and Z { c , do not depend on c, and fiiidm = c eq + L s . 
On the other hand, obviously, c cq e f^idm- Therefore, procedures of solving 
of Eq. (27) can be initiated by choosing the linear approximation, f2^ m = 
c eq + L s C c q , in the neighborhood of the equilibrium c eq , and then continued 
parametrically into the nonlinear domain. Computational technologies of a 
continuation of solutions with respect to parameters are well developed (see, 
for example, Khibnik, Kuznetsov, Levitin fc Nikolaev (1993 ); [Roose fc Spencc 



(1990|) ). The problem of the relevant parameterization is solved locally: In the 
neighborhood of a given point c° one can choose Z s c (c — c°) for a characteri- 
zation of the vector c. In this case, the space of parameters is L s c . There exist 
other, physically motivated ways to parameterize manifolds ( porban fc Karliri 
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(1992); see also section 4.1 below). 



There are two drawbacks of the ILDM method which call for its refinement: 
(i) "Intrinsic" does not imply "invariant" . Eq. (27) is not invariant of the 
dynamics (4). If one differentiates Eq. (27) in time due to Eq. (4), one obtains 
a new equation which is the implication of Eq. (27) only for linear systems. In 
a general case, the motion c(t) takes off the findm- Invariance of a manifold 
Q means that «7(c) touches f2 in every point c e S7. It remains unclear how 
the ILDM (27) corresponds with this condition. Thus, from the dynamical 
perspective, the status of the ILDM remains not well defined, or "ILDM is 
ILDM", defined self-consistently by Eq. (27), and that is all what can be 
said about it. (ii) From the geometrical standpoint, spectral decomposition of 
Jacobian fields is not the most attractive way to compute manifolds. If we are 
interested in the behavior of trajectories, how they converge or diverge, then 
one should consider the symmetrized part of F c , rather than F c itself. 

Symmetric part, F s ^ m = (1/2)(Fq + F c ), defines the dynamics of the dis- 
tance between two solutions, c and c', in a given local Euclidean metrics. 
Skew-symmetric part defines rotations. If we want to study manifolds based 
on the argument about convergence/divergence of trajectories, then we should 
use in Eq. (27) the spectral projector Z f ^ m for the operator Fq 111 . This, by 
the way, is also a significant simplification from the standpoint of computa- 
tions. It remains to choose the metrics. This choice is unambiguous from the 
thermodynamic perspective. In fact, there is only one choice which fits into 
the physical meaning of the problem, this is the metrics associated with the 
thermodynamic (or entropic) scalar product, 



where He is the matrix of second-order derivatives of G (7). In the equilib- 
rium, operator F c ^ is selfajoint with respect to this scalar product (Onsager's 
reciprocity relations). Therefore, the behavior of the ILDM in the vicinity of 
the equilibrium does not alter under the replacement, F c ^ = F s ^. In terms 
of usual matrix representation, we have: 



where F c is the ordinary transposition. 

The ILDM constructed with the help of the symmetrized Jacobian field will be 
termed the symmetric entropic intrinsic low-dimensional manifold (SEILDM). 
Selfadjointness of F s ^ m (29) with respect to the thermodynamic scalar product 
(28) simplifies considerably computations of spectral decomposition. More- 
over, it becomes necessary to do spectral decomposition in only one point - 



{{x,y)) = (x,H c y), 



(28) 



F s r^l(F c + H^F T c H c ), 



(29) 
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in the equilibrium. Perturbation theory for selfadjoint operators is a very well 
developed subject ( [Kato, 19761) , which makes it possible to easily extend the 
spectral decomposition with respect to parameters. A more detailed discussion 
of the selfadjoint linearization will be given below in section 5.2. 

Thus, when the geometric picture behind the decomposition of motions is 
specified, the physical significance of the ILDM becomes more transparent, 
and it leads to its modification into the SEILDM. This also gains simplicity 
in the implementation by switching from non- selfadjoint spectral problems to 
selfadjoint. The quantitative estimate of this simplification is readily available: 
Let d be the dimension of the phase space, and k the dimension of the ILDM 
(k = dimL s c ). The space of all the projectors Z with the fc-dimensional image 
has the dimension D = 2k(d — k). The space of all the selfadjoint projectors 
with the fc-dimensional image has the dimension D sym = k(d — k). For d = 20 
and k = 3, we have D = 102 and _D sym = 51. When the spectral decomposition 
by means of parametric extension is addressed, one considers equations of the 
form: 

^ Few, VFcm), (30) 



where r is the parameter, and V-Fc = W«7(c) is the differential of the 
Jacobian field. For the selfadjoint case, where we use = F s c ym instead of F c , 
this system of equations has twice less independent variables, and also the 
right hand is of a simpler structure. 

It is more difficult to improve on the first of the remarks (ILDM is not invari- 
ant). The following naive approach may seem possible: 

(i) Take fludm = cCq + L s C c q in a neighborhood U of the equilibrium c eq . [This 
is also a useful initial approximation for solving Eq. (27)]. 

(ii) Instead of computing the solution to Eq. (27), integrate the kinetic equa- 
tions (4) backwards in the time. It is sufficient to take initial conditions c(0) 
from a dense set on the boundary, dU D (c eq + L|,oq), and to compute solutions 
c(t), t < 0. 

(iii) Consider the obtained set of trajectories as an approximation of the slow 
invariant manifold. 

This approach will guarantee invariance, by construction, but it is prone to 
pitfalls in what concerns the slowness. Indeed, the integration backwards in the 
time will see exponentially divergent trajectories, if they were exponentially 
converging in the normal time progress. This way one finds some invariant 
manifold which touches c cq + L 8 C eq in the equilibrium. Unfortunately, there 
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are infinitely many such manifolds, and they fill out almost all the space of 
concentrations. However, we must select the slow component of motions. Such 
a regularization is possible. Indeed, let us replace in Eq. (4) the vector field 
•7(c) by the vector field Z s ^ ym J(c), and obtain a regularized kinetic equation, 

c = Z s * ym J(c). (31) 



Let us replace integration backwards in time of the kinetic equation (4) in the 
naive approach described above by integration backwards in time of the regu- 
larized kinetic equation (31). With this, we obtain a rather convincing version 
of the ILDM (SEILDM). Using Eq. (30), one also can write down an equation 
for the projector Z s ^ ym , putting r = t. Replacement of Eq. (4) by Eq. (31) also 
makes the integration backwards in time in the naive approach more stable. 
However, regularization will again conflict with invariance. The "naive refine- 
ment" after the regularization (31) produces just a slightly different version of 
the ILDM (or SEILDM) but it does not construct the slow invariant manifold. 
So, where is the way out? We believe that the ILDM and its version SEILDM 
are, in general, good initial approximations of the slow manifold. However, if 
one is indeed interested in finding the invariant manifold, one has to write out 
the true condition of invariance and solve it. As for the initial approximation 
for the method of invariant manifold one can use any ansatz, in particular, 
the SEILDM. 

The problem of a complete decomposition of kinetic equations can be solved 
indeed in some cases. The first such solution was the spectral decomposition 
for linear systems ( |Wei fc Prater, 1962| ). Decomposition is sometimes possible 
also for nonlinear systems ([Li, Rabitz fc Toth (1994] ); |T6th, Li, Rabitz fcj 
Tomlin (19971 )). The most famous example of a complete decomposition of 
infinite-dimensional kinetic equation is the complete integrability of the space- 
independent Boltzmann equation for Maxwell's molecules found by |Bobylev 



T988|). However, in a general case, there exist no analytical, not even a twice 



differentiable transformation which would decouple modes. The well known 
Grobman-Hartman theorem ( jHartman, 1963| , |1982| ) states only the existence 



of a continuous transform which decomposes modes in a neighborhood of the 
equilibrium. For example, the analytic planar system, dx/dt = —x, dy/dt = 
—2y + x 2 , is not C 2 linearizable. These problems remain of interest (|Chicone,| 
200(J| ). Therefore, in particular, it becomes quite ineffective to construct such 
a transformation in a form of a series. It is more effective to solve a simpler 
problem of extraction of a slow invariant manifold QBeyn fc Kless, 1998| ). 



Sensitivity analysis ( [Rabitz, Kramer fc Dacol, 1983 ; Rabitz, 1987 ; Lam & 



Goussis, 1994| ) makes it possible to select essential variables and reactions, 
and to decompose motions into fast and slow. In a sense, the ILDM method 
is a development of the sensitivity analysis. Recently, a further step in this 
direction was done by |Zhu fc Petzold (1999| ). In this work, the authors use 
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a nonlocal in time criterion of closeness of solutions of the full and of the 
reduced systems of chemical kinetics. They require not just a closeness of 
derivatives but a true closeness of the dynamics. 

Let us be interested in the dynamics of the concentrations of just a few species, 
Ai, . . . , A p , whereas the rest of the species, A p+1 , . . . , A n are used for build- 
ing the kinetic equation, and for understanding the process. Let c goa i be the 
concentration vector with components Ci,...,c p , c goa i(t) be the correspond- 
ing components of the solution to Eq. (4), and be the solution to the 
simplified model with corresponding initial conditions. Zhu fc Petzold (1999| ) 



suggest to minimize the difference between c goa i(t) and on the segment 
t G [0, T]: ||c goa i(i) — Cg^H — > min. In the course of the optimization under 
certain restrictions one selects the optimal (or appropriate) reduced model. 
The sequential quadratic programming method and heuristic rules of sorting 
the reactions, substances etc were used. In the result, for some stiff systems 
studied, one avoids typical pitfalls of the local sensitivity analysis. In simpler 
situations this method should give similar results as the local methods. 



2. 7 Thermodynamic criteria for selection of important reactions 



One of the problems addressed by the sensitivity analysis is the selection of 
the important and discarding the unimportant reactions. |Bykov, Yablonskii 



fc Akramov (1977|) suggested a simple principle to compare importance of dif- 



ferent reactions according to their contribution to the entropy production (or, 
which is the same, according to their contribution to G). Based on this princi- 
ple, |Dimitrov (1982|) described domains of parameters in which the reaction of 



hydrogen oxidation, H2 + O2 + M, proceeds due to different mechanisms. For 
each elementary reaction, he has derived the domain inside which the contri- 
bution of this reaction is essential (nonnegligible) . Due to its simplicity, this 
entropy production principle is especially well suited for analysis of complex 
problems. In particular, recently, a version of the entropy production princi- 
ple was used in the problem of selection of boundary conditions for Grad's 
moment equations flStruchtrup fc Weiss, 199~B| |Grmela, Karlin fc Zmievsk" 



2002j ). For ideal systems (9), the contribution of the sth reaction to G has a 



particularly simple form: 



(W + \ ■ J-^ ■ 
G s = -W s \n(^-\, G = Y,G S . (32) 



For nonideal systems, the corresponding expressions (13) are also not too 
complicated. 
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3 Outline of the method of invariant manifold 

In many cases, dynamics of the d- dimensional system (4) leads to a manifold of 
a lower dimension. Intuitively, a typical phase trajectory behaves as follows: 
Given the initial state c(0) at t — 0, and after some period of time, the 
trajectory comes close to some low- dimensional manifold Q, and after that 
proceeds towards the equilibrium essentially along this manifold. The goal is 
to construct this manifold. 

The starting point of our approach is based on a formulation of the two main 
requirements: 

(i) . Dynamic invariance: The manifold fl should be (positively) invariant un- 
der the dynamics of the originating system (4): If c(0) G f2, then c(t) G ft for 
each t > 0. 

(ii) . Thermodynamic consistency of the reduced dynamics: Let some (not oblig- 
atory invariant) manifold Q is considered as a manifold of reduced description. 
We should define a set of linear operators, Pc, labeled by the states c G f2, 
which project the vectors «7(c), c G Q onto the tangent bundle of the manifold 
f2, thereby generating the induced vector field, PcJ(c), c G fl. This induced 
vector field on the tangent bundle of the manifold fl is identified with the 
reduced dynamics along the manifold f2. The thermodynamicity requirement 
for this induced vector field reads 



In order to meet these requirements, the method of invariant manifold suggests 
two complementary procedures: 

(i) . To treat the condition of dynamic invariance as an equation, and to solve 
it iteratively by a Newton method. This procedure is geometric in its nature, 
and it does not use the time dependence and small parameters. 

(ii) . Given an approximate manifold of reduced description, to construct the 
projector satisfying the condition (33) in a way which does not depend on the 
vector field J. 

We shall now outline both these procedures starting with the second. The so- 
lution consists, in the first place, in formulating the thermodynamic condition 
which should be met by the projectors Pc'. For each c G fl, let us consider 
the linear functional 




(33) 




(34) 
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Then the thermodynamic condition for the projectors reads: 

keri-*c Q kerM£, for each c G f2. 



(35) 



Here kerPc is the null space of the projector, and kerM£ is the hyperplane 
orthogonal to the vector M£. It has been shown flGorban fc Karlin, 1992| , [L994| ) 



that the condition (35) is the necessary and sufficient condition to establish 
the thermodynamic induce vector field on the given manifold f2 for all possible 
dissipative vector fields J simultaneously 

Let us now turn to the requirement of invariance. By a definition, the manifold 
f2 is invariant with respect to the vector field J if and only if the following 
equality is true: 

[1 - P] J(c) = 0, for each c E ft. (36) 



In this expression P is an arbitrary projector on the tangent bundle of the 
manifold f2. It has been suggested to consider the condition (36) as an equation 
to be solved iteratively starting with some appropriate initial manifold. 

Iterations for the invariance equation (36) are considered in the section 5. The 
next section presents construction of the thermodynamic projector using a 
specific parameterization of manifolds. 



4 Thermodynamic projector 



4-1 Thermodynamic parameterization 



In this section, f2 denotes a generic p-dimensional manifold. First, it should 
be mentioned that any parameterization of Q generates a certain projector, 
and thereby a certain reduced dynamics. Indeed, let us consider a set of m 
independent functionals M(c) = {Mi(c), . . . , M p (c)}, and let us assume that 
they form a coordinate system on f2 in such a way that f2 = c(M), where 
c(M) is a vector function of the parameters Mi, . . . , M p . Then the projector 
associated with this parameterization reads: 

p dc(M) 

p c(M )X = £ -^T7 ± Nr j \M){VM j \ c{M) ,x), (37) 

i,j=l % 
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where N^ 1 is the inverse to the p x p matrix: 

N(M) = WiVM^dc/dMj^. (38) 

This somewhat involved notation is intended to stress that the projector (37) 
is dictated by the choice of the parameterization. Subsequently, the induced 
vector field of the reduced dynamics is found by applying projectors (37) on 
the vectors J(c(M)), thereby inducing the reduced dynamics in terms of the 
parameters M as follows: 

M t = ]TA^(M)(VM, \ c{M) ,J(c(M))), (39) 



Depending on the choice of the parameterization, dynamic equations (39) 
are (or are not) consistent with the thermodynamic requirement (33). The 
thermodynamic parameterization makes use of the condition (35) in order to 
establish the thermodynamic projector. Specializing to the case (37), let us 
consider the linear functionals, 

DM t \ c{M) (x) = (VM t \ c{M) ,x}. (40) 



Then the condition (35) takes the form: 

p 

r 

i=i 



ftkevDM, | cw CkerM* (M) , (41) 



that is, the intersection of null spaces of the functionals (40) should belong to 
the null space of the differential of the Lyapunov function G, in each point of 
the manifold f2. 

In practice, in order to construct the thermodynamic parameterization, we 
take the following set of functionals in each point c of the manifold f2: 



M l (x) = M* c (x), cefl (42) 
M i (x) = (m i ,x), i = 2,...,p (43) 

It is required that vectors VG(c), m 2 , . . . , m p are linearly independent in each 
state c e f2. Inclusion of the functionals (34) as a part of the system (42) and 
(43) implies the thermodynamic condition (41). Also, any linear combination 
of the parameter set (42), (43) will meet the thermodynamicity requirement. 

It is important to notice here that the thermodynamic condition is satisfied 
whatsoever the functionals M 2 , . . . , M p are. This is very convenient for it gives 
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an opportunity to take into account the conserved quantities correctly. The 
manifolds we are going to deal with should be consistent with the conservation 
laws (2). While the explicit characterization of the phase space V is a problem 
on its own, in practice, it is customary to work in the n-dimensional space 
while keeping the constraints (2) explicitly on each step of the construction. 
For this technical reason, it is convenient to consider manifolds of the dimen- 
sion p > I, where I is the number of conservation laws, in the n-dimensional 
space rather than in the phase space V. The thermodynamic parameterization 
is then concordant also with the conservation laws if I of the linear functionals 
(43) are identified with the conservation laws. In the sequel, only projectors 
consistent with conservation laws are considered. 

Very frequently, the manifold f2 is represented as a p-parametric family c(a\, . . . ,a p ), 
where a, are coordinates on the manifold. The thermodynamic re-parameterization 
suggests a representation of the coordinates in terms of M£, M 2 , . . . , M p 
(42), (43). While the explicit construction of these functions may be a formidable 
task, we notice that the construction of the thermodynamic projector of the 
form (37) and of the dynamic equations (39) is relatively easy because only 
the derivatives dc/dMi enter these expressions. This point was discussed in a 
detail by |Gorban fc Karlin (1992j [199g) . 



4-2 Decomposition of motions: Thermodynamics 



Finally, let us discuss how the thermodynamic projector is related to the de- 
composition of motions. Assuming that the decomposition of motions near 
the manifold fi is true indeed, let us consider states which were initially close 
enough to the manifold f2. Even without knowing the details about the evo- 
lution of the states towards f2, we know that the Lyapunov function G was 
decreasing in the course of this evolution. Let us consider a set of states U c 
which contains all those vectors c! that have arrived (in other words, have been 
projected) into the point c 6 f2. Then we observe that the state c furnishes 
the minimum of the function G on the set Uc- If a state c! e Uc, and if it 
deviates small enough from the state c so that the linear approximation is 
valid, then c' belongs to the affine hyperplane 

T c = c + ker M* , cefl. (44) 



This hyperplane actually participates in the condition (35). The consideration 
was entitled 'thermodynamic' ( |Gorban fc Karlin, 1992| ) because it describes 
the states c 6 f2 as points of minimum of the function G over the correspond- 
ing hyperplanes (44). 
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5 Corrections 



5.1 Preliminary discussion 

The thermodynamic projector is needed to induce the dynamics on a given 
manifold in such a way that the dissipation inequality (33) holds. Coming back 
to the issue of constructing corrections, we should stress that the projector 
participating in the invariance condition (36) is arbitrary. It is convenient to 
make use of this point: When Eq. (36) is solved iteratively, the projector may 
be kept non-thermodynamic unless the induced dynamics is explicitly needed. 

Let us assume that we have chosen the initial manifold, Qq, together with 
the associated projector Pq, as the first approximation to the desired man- 
ifold of reduced description. Though the choice of the initial approximation 
fl depends on the specific problem, it is often reasonable to consider quasi- 
equilibrium or quasi steady-state approximations. In most cases, the manifold 
S7o is n °t an invariant manifold. This means that Qq does not satisfy the 
invariance condition (36): 

A = [1 - P ] J(c ) ^ 0, for some c G Q . (45) 

Therefore, we seek a correction c\ = c + 5c. Substituting P = P and 
c = c + 5c into the invariance equation (36), and after the linearization in 
5c, we derive the following linear equation: 

[1 - P ] [J(c) + LcM = 0, (46) 

where L Co is the matrix of first derivatives of the vector function «/, computed 
in the state c e fio- The system of linear algebraic equations (46) should be 
supplied with the additional condition. 

P 5c = 0. (47) 

In order to illustrate the nature of the Eq. (46), let us consider the case of 
linear manifolds for linear systems. Let a linear evolution equation is given 
in the finite-dimensional real space: c = Lc, where L is negatively definite 
symmetric matrix with a simple spectrum. Let us further assume the quadratic 
Lyapunov function, G(c) = (c, c). The manifolds we consider are lines, 1(a) = 
ae, where e is the unit vector, and a is a scalar. The invariance equation for 
such manifolds reads: e(e, Le)—Le = 0, and is simply a form of the eigenvalue 
problem for the operator L. Solutions to the latter equation are eigenvectors 
ej, corresponding to eigenvalues Aj. 
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Assume that we have chosen a line, Iq = aeo, defined by the unit vector eo, and 
that eo is not an eigenvector of L. We seek another line, li = ae±, where e x is 
another unit vector, ei = yi/\\yi\\, Hi = e^+Sy. The additional condition (47) 
now reads: (5y, e ) = 0. Then the Eq. (46) becomes [1 — e (e , -)]L[e +5y] = 0. 
Subject to the additional condition, the unique solution is as follows: e + 5y = 
(eo, L~ 1 eo)~ 1 L~ 1 eo. Rewriting the latter expression in the eigen-basis of L, 
we have: eo + 5y oc J2i ^i le i{ e i-> e o)- The leading term in this sum corresponds 
to the eigenvalue with the minimal absolute value. The example indicates that 
the method of linearization (46) seeks the direction of the slowest relaxation. 
For this reason, the method (46) can be recognized as the basis of an iterative 
method for constructing the manifolds of slow motions. 

For the nonlinear systems, the matrix L Co in the Eq. (46) depends nontrivially 
on c . In this case the system (46) requires a further specification which will 
be done now. 



5.2 Symmetric linearization 

The invariance condition (36) supports a lot of invariant manifolds, and not all 
of them are relevant to the reduced description (for example, any individual 
trajectory is itself an invariant manifold). This should be carefully taken into 
account when deriving a relevant equation for the correction in the states of 
the initial manifold f2 which are located far from equilibrium. This point 
concerns the procedure of the linearization of the vector field J, appearing 
in the equation (46). We shall return to the explicit form of the Marcelin-Dc 
Donder kinetic function (11). Let c is an arbitrary fixed element of the phase 
space. The linearization of the vector function J (12) about c may be written 
J(c + 5c) PS J(c) + L c 5c where the linear operator L c acts as follows: 

L c x = ir~f s [Ws + (c)(<x s ,H c x) - W s -(c)((3 s ,H c x)}. (48) 

s=l 

Here H c is the matrix of second derivatives of the function G in the state c 
[see Eq. (7)]. The matrix L c in the Eq. (48) can be decomposed as follows: 

L c = L' c + L" c . (49) 
Matrices L' c and L" c act as follows: 

L 'c* = ~\ E[^ s + (c) + W7(c)] 7 .<7., Hex), (50) 

Z s=l 
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L'ix = \ £[W^(c) - W;{c)] ls {cx s + fi„ H c x). (51) 

Z 8=1 

Some features of this decomposition are best seen when we use the thermo- 
dynamic scalar product (28): The following properties of the matrix L' c are 
verified immediately: 

(i) The matrix L' c is symmetric in the scalar product (28): 

((x,L' c y)) = ((y,L' c x)). (52) 

(ii) The matrix L' c is nonpositive definite in the scalar product (28): 

{(x,L' c x))<0. (53) 

(iii) The null space of the matrix L' c is the linear envelope of the vectors 
H^bi representing the complete system of conservation laws: 

keiL' c = LmiH^bi, i = 1, . . . , 1} (54) 

(iv) If c = c cc i, then W+(c eq ) = W7(c eq ), and 

Thus, the decomposition Eq. (49) splits the matrix Lc in two parts: one 
part, Eq. (50) is symmetric and nonpositive definite, while the other part, 
Eq. (51), vanishes in the equilibrium. The decomposition Eq. (49) explicitly 
takes into account the Marcelin-De Donder form of the kinetic function. For 
other dissipative systems, the decomposition (49) is possible as soon as the 
relevant kinetic operator is written in a gain-loss form [for instance, this is 
straightforward for the Boltzmann collision operator]. 

In the sequel, we shall make use of the properties of the operator L' c (50) for 
constructing the dynamic correction by extending the picture of the decom- 
position of motions. 

5. 3 Decomposition of motions: Kinetics 

The assumption about the existence of the decomposition of motions near the 
manifold of reduced description f2 has led to the thermodynamic specifica- 
tions of the states c e f2. This was accomplished in the section 4.2, where the 
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(55) 



thermodynamic projector was backed by an appropriate variational formula- 
tion, and this helped us to establish the induced dynamics consistent with 
the dissipation property. Another important feature of the decomposition of 
motions is that the states c G f2 can be specified kinetically. Indeed, let us do 
it again as if the decomposition of motions were valid in the neighborhood of 
the manifold Q, and let us 'freeze' the slow dynamics along the S7, focusing on 
the fast process of relaxation towards a state c G f2. From the thermodynamic 
perspective, fast motions take place on the afline hyperplane c + 5c G Tq , 
where r Co is given by Eq. (44). From the kinetic perspective, fast motions on 
this hyperplane should be treated as a relaxation equation, equipped with the 
quadratic Lyapunov function 5G = {{5c, 5c)), Futhermore, we require that the 
linear operator of this evolution equation should respect Onsager's symme- 
try requirements (selfadjointness with respect to the entropic scalar product). 
This latter crucial requirement describes fast motions under the frozen slow 
evolution in the similar way, as all the motions near the equilibrium. 

Let us consider now the manifold Qq which is not the invariant manifold of the 
reduced description but, by our assumption, is located close to it. Consider a 
state c G fl , and the states c + 5c close to it. Further, let us consider an 
equation 



Due to the properties of the operator L' c (50), this equation can be regarded 
as a model of the assumed true relaxation equation near the true manifold of 
the reduced description. For this reason, we shall use the symmetric operator 
L' c (50) instead of the linear operator L c when constructing the corrections. 

5.4 Symmetric iteration 

Let the manifold Qq and the corresponding projector Pq are the initial ap- 
proximation to the invariant manifold of the reduced description. The dynamic 
correction c± = c + 5c is found upon solving the following system of linear 
algebraic equations: 



Here L' Cq is the matrix (50) taken in the states on the manifold Qq. An im- 
portant technical point here is that the linear system (57) always has the 
unique solution for any choice of the manifold Q. This point is crucial since it 
guarantees the opportunity of carrying out the correction process for arbitrary 
number of steps. 



5 c = L' n 5c. 



(56) 



[1 - P ] J (cq) + L' Co 5c = 0, P 5c = 0. 



(57) 
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6 The method of invariant manifold 



We shall now combine together the two procedures discussed above. The re- 
sulting method of invariant manifold intends to seek iteratively the reduced 
description, starting with an initial approximation. 

(i). Initialization. In order to start the procedure, it is required to choose the 
initial manifold Qq, and to derive corresponding thermodynamic projector 
Pq. In the majority of cases, initial manifolds are available in two different 
ways. The first case are the quasi-equilibrium manifolds described in the sec- 
tion 2.3. The macroscopic parameters are Mj = q = (rrij, c), where rrij is 
the unit vector corresponding to the specie A+. The quasi-equilibrium mani- 
fold, Co (Mi, . . . , Mfc, -E>i, . . . , Bi), compatible with the conservation laws, is the 
solution to the variational problem: 



G — > min , (mj, c) = q, i = 1, . . . , k, (58) 
(bj,c) = Bj, j = l,...,l. 

In the case of quasi-equilibrium approximation, the corresponding thermody- 
namic projector can be written most straightforwardly in terms of the variables 

Pox = £ ^(m, x) + £ x). (59) 

i=l 1 i=l 1 



For quasi-equilibrium manifolds, a reparameterization with the set (42), (43) 
is not necessary QGorban fc Karlin (1992| ); |Gorban fc Karlin (1994|) ). 



The second source of initial approximations are quasi-stationary manifolds 
(section 2.5). Unlike the quasi-equilibrium case, the quasi-stationary manifolds 
must be reparameterized in order to construct the thermodynamic projector. 

(ii). Corrections. Iterations are organized in accord with the rule: If c m is 
the mth approximation to the invariant manifold, then the correction c m+1 = 
c m + Sc is found from the linear algebraic equations, 



[1 - P m }(J(c m ) + L' c Jc) = 0, (60) 
P m Sc = 0. (61) 

Here L' Cm is the symmetric matrix (50) evaluated at the mth approximation. 
The projector P m is not obligatory thermodynamic at that step, and it is 
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taken as follows: 

k dc m 1 dc 

P m x = Y,^( m ^ x ) + Y,-Q^( b ^ x )- ( 62 ) 

i=l 1 i=l 1 

(iii). Dynamics. Dynamics on the mth manifold is obtained with the thermo- 
dynamic re-parameterization. 

In the next section we shall illustrate how this all works. 
7 Illustration: Two-step catalytic reaction 

Here we consider a two-step four-component reaction with one catalyst A 2 : 

A Y + A 2 ^ A 3 ^A 2 + A 4 . (63) 

We assume the Lyapunov function of the form (9), G — J2t=i Cj[l n (ci/c° q ) — 1]. 
The kinetic equation for the four-component vector of concentrations, c = 
(ci, c 2 , c 3 , c 4 ), has the form 

c = ll W l + l2 W 2 . (64) 

Here ~f l 2 are stoichiometric vectors, 

7l = (-1,-1, 1,0), 72 = (0,1, -1,1), (65) 

while functions W^ 2 are reaction rates: 

W 1 = k\c x c 2 - fcf c 3 , 1^2 = A;Jc 3 - k 2 c 2 c 4 . (66) 

Here /c^ 2 are reaction rate constants. The system under consideration has two 
conservation laws, 

ci + c 3 + c 4 = Si, c 2 + c 3 = 5 2 , (67) 

or (6 lj2 ,c) = B 1)2 , where &i = (1,0,1,1) and b 2 = (0,1,1,0). The nonlinear 
system (64) is effectively two-dimensional, and we consider a one-dimensional 
reduced description. 

We have chosen the concentration of the specie A 1 as the variable of reduced 
description: M = Ci, and c\ = (m, c), where m = (1,0,0,0). The initial 
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manifold Cq(M) was taken as the quasi-equilibrium approximation, i.e. the 
vector function c is the solution to the problem: 

G — > min for (m, c) — c±, (b±, c) = B\, (b 2 , c) = B 2 . (68) 



The solution to the problem (68) reads: 



c 01 = d, (69) 
c 2 = B 2 - 0(ci), 

C O 3 = 0(Cl), 

004 = ^1 - Ci - 0(Ci), 

0(M) = A( Cl ) - y /A 2 (c 1 )-B 2 (B 1 -c 1 ), 

A( ^ B^-c^ + cfjcT + cf-c,) 
A{ Cl ) - ^ . 

The thermodynamic projector associated with the manifold (69) reads: 

P x = ^(m, x) + + ||(6 2 , x). (70) 



Computing A = [1 — P ]«^( c o) we find that the inequality (45) takes place, 
and thus the manifold Cq is not invariant. The first correction, C\ = c + 5c, 
is found from the linear algebraic system (60) 



(l-P )L' 8c=-[l-Po]J(co), (71) 

5d = 
5ci + 5c3 + 5c4 = 

5c 3 + 5c 2 = 0, (72) 

where the symmetric 4x4 matrix L' has the form (we write instead of c 
in the subscript in order to simplify notations): 

W+(co) + Wr(co) 7ii W 2 + (c ) + W 2 (c ) l2l 
2 c / 2 cqi 



The explicit solution Ci(ci, B x , B 2 ) to the linear system (71) is easily found, 
and we do not reproduce it here. The process was iterated. On the k + 1 
iteration, the following projector P k was used: 
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Notice that projector P k (74) is the thermodynamic projector only if k = 0. As 
we have already mentioned it above, in the process of finding the corrections 
to the manifold, the non-thermodynamic projectors are allowed. The linear 
equation at the k + 1 iteration is thus obtained by replacing c , Pq, and L' Q 
with Cfc, Pfc, and L' k in all the entries of the Eqs. (71) and (73). 

Once the manifold c& was obtained on the fcth iteration, we derived the corre- 
sponding dynamics by introducing the thermodynamic parameterization (and 
the corresponding thermodynamic projector) with the help of the function 
(42). The resulting dynamic equation for the variable C\ in the fcth approxi- 
mation has the form: 

(VG l^dc/dcjc! = (VG | Cfc , J(c k )). (75) 
Here [VG | C J, = Hc ki /c?]. 

Analytic results were compared with the results of the numerical integration. 
The following set of parameters was used: 

kf = 1.0, fcf = 0.5, k+ = 0.4, = 1.0; 
c f = 0.5, c c 2 q = 0.1, c c 3 q = 0.1, = 0.4, 
B x = 1.0, B 2 = 0.2. 

Direct numerical integration of the system has demonstrated that the manifold 
C3 = in the plane (01,03) attracts all individual trajectories. Thus, the 
reduced description in this example should extract this manifold. 

Fig. 1 demonstrates the quasi-equilibrium manifold (69) and first two correc- 
tions found analytically. It is apparent that while the initial quasi-equilibrium 
approximation is in a poor agreement with the reduced description, the cor- 
rections rapidly improve the situation. This confirms our expectation of an 
advantage of using iteration methods in comparison to methods based on a 
small parameter expansions. 



8 Method of invariant manifold without a priori parameterization 

Formally, the method of invariant manifold does not require a global parame- 
terization of the manifolds. However, in most of the cases, one makes use of a 
priori defined "macroscopic" variables M. This is motivated by the choice of 
quasi-equilibrium initial approximations. 

Let a manifold fl be defined in the phase space of the system, its tangent space 
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in the point c be Tcfi. How to define the projector of the whole concentrations 
space onto T c f2 without using any a priori parameterization of ft? 

The basis of the answer to this question is the condition of thermodynamicity 
(35). Let us denote E as the concentration space, and consider the problem 
of the choice of the projector in the quadratic approximation to the thermo- 
dynamic potential G: 

G q = (g, H c Ac) + i(Ac, H c Ac) = ((g, Ac}} + I ((Ac, Ac)), (76) 

where He is the matrix of the second-order derivatives of G (7), g — H^'VG, 
Ac is the deviation of the concentration vector from the expansion point. 

Let a linear subspace T be given in the concentrations space E. Problem: For 
every Ac + T, and for every g G E, define a subspace Lac such that: (i) Lac 
is a complement of T in E: 

L AC + T = E, L Ac nT = {0}. 



(ii) Ac is the point of minimum of G q on Lac + Ac: 



Ac = arg^mm G^x). (77) 



Besides (i) and (ii), we also impose the requirement of a maximal smoothness 
(analyticity) on Lac as a function of g and Ac. Requirement (77) implies that 
Ac is the quasi-equilibrium point for the given Lac, while the problem in a 
whole is the inverse quasi-equilibrium problem: We construct Lac such that 
T will be the quasi-equilibrium manifold. Then subspaces Lac will actually 
be the kernels of the quasi-equilibrium projector. 

Let f 1 , . . . , f k be the orthonormalized with respect to ((•, •)) scalar product 
basis of T, vector h be orthogonal to T, ((h, h}} — 1, g — af 1 + f3h. Condition 
(77) implies that the vector VG is orthogonal to Lac m the point Ac. 

Let us first consider the case (3 = 0. The requirement of analyticity of Lac 
as the function of a and Ac implies Lac — Lq + o(l), where Lq = T L is 
the orthogonal completement of T with respect to scalar product ((•,•)). The 
constant solution, Lac = Lq also satisfies (77). Let us fix a ^ 0, and extend 
this latter solution to (3 ^ 0. With this, we obtain a basis, l± : . . . , Z„_ fc . Here 
is the simplest construction of this basis: 

(3f 1 -ja + A Cl )h 
1 (/? 2 + (a + Ac 1 ) 2 ) 1 /2' ^ 
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where Aci = ((Ac, is the first component in the expansion, Ac = 
X^Acj/j. The rest of the basis elements, I2, . . . ,l n -k form the orthogonal 
completement of T© (h) with respect to scalar product ((•, •)), (h) is the line 
spanned by h. 

Dependence Lac (78) on Ac, a and f3 is singular: At a + Aci, vector Zi G 
T, and then Lac is not the completement of T in L anymore. For a / 0, 
dependence Lac gives one of the solutions to the inverse quasi-equilibrium 
problem in the neighborhood of zero in T. We are interested only in the limit, 



{Pfi - n/i 



lim L AC = LhW T == , h, ■ ■ • , h-k \- (79) 
ac^O [ V" + P J 



Finally, let us define now the projector Pc of the space E onto Tc^l- If 
i3" c 1 VG G T c f2, then P c is the orthogonal projector with respect to the 
scalar product ((•,•)): 

k 

Pcz^MUz)). (80) 
i=i 



If H^VG <£ T c n, then, according to Eq. (79), 



where . . . , is the orthonormal with respect to ((•, •)) basis of T c f2, 
h is orthogonal to T, = 1, H^VG = af 1 + Z x = ((3f 1 - 

ah)/^+W, {(f l M)) = p/V^ r TW- 

Thus, for solving the invariance equation iteratively, one needs only projector 
Pc (81), and one does not need a priori parameterization of f2 anymore. 



9 Method of invariant grids 



Grid-based approximations of manifolds are attractive from the computational 
perspective. Since no a priori parameterization is required in the method of 
invariant manifold, in this section we develop its grid-based realization. Let 
us consider a regular grid Q in R k , and its mapping F into the concentrations 
space E. It makes sense to consider only F which map a finite part of the grid 
into the phase space V. This part of the map is termed essential. Extension of 
the map F onto the rest of the nodes is done by a simple (for example, linear) 
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extrapolation of the essential part (in practice, one needs to extrapolate only 
onto the next neighbors of the essential nodes). 

Let operators of grid differentiation Di where be defined for functions on the 
grid, where i = 1, . . . , k label grid coordinates Xj. With this, the tangent space 
to the image of the grid in the point c(x) = F(x) is defined for each node of 
the grid x: 



T x = Lin{</?i, . . .,(pk}, 
(Pi = Dic(x) = (Didix), . . . , DiC n (x)). (82) 

The grid is termed invariant if, for each essential node, 

J{c{x)) e T x . 



For the essential nodes, we write down the invariance equation with the pro- 
jector, Pc{x) '■ E — > T x : This equation is solved using the Newton method 
as it was described above in the section 6. A good initial approximation is a 
linear map of the grid on the affine manifold corresponding to slow relaxation 
in the vicinity of the equilibrium. It is convenient to take this map isometric 
with respect to the metrics generated by the entropic scalar product in the 
equilibrium. 

If the vector field of the reduced model, c = P C ( x )J(c(x)), is defined on the 
nodes F(x), then one can define the dynamics X{ on the nodes. In order to 
do this, we expand c over ipf c = Yjl=i a>ifi- The dynamics on the nodes is 
then defined by equations, Xi = di. Using interpolation, we can define the 
vector field x within the essential cells of the grid (those cells for which the all 
the nodes are essential). The system of equations thus obtained models the 
dynamics on the invariant manifold. 

The essence of this construction is that, by solving a set of uncomplicated 
linear equations arising from linearization of the invariance equations on the 
nodes one gets a reliable numerical scheme for constructing invariant mani- 
folds. The use of the grid differentiation rather than a different iable approxi- 
mation to the manifold makes the scheme suited for parallel realizations. We 
stress it once again that such realizations are only possible if no a priori global 
parameterization of manifolds is required. Further refinements of the scheme, 
taking into account the process of moving the inessential nodes into the phase 
space, and the opposite process of essential nodes leaving the phase space can 
be done based in the same way as for grid-based data analysis ( |Gorban~^ 
Kossiev, 1999|; |Gorban fc Zinovyev, 2(J(JT|). 
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10 Method of invariant manifold for open systems 



One of the problems to be focused on when studying closed systems is to 
prepare extensions of the result for open or driven by flows systems. External 
flows are usually taken into account by addidinal terms in the kinetic equations 
(4): 

c = J( c ) + n. (83) 



Zero- order approximation assumes that the flow does not change the invariant 
manifold. Equations of the reduced dynamics, however, do change: Instead of 
J(c(M)) we substitute J(c(M)) + II into Eq. (39): 

Mi = t^ 1 ^!^, J(c(M)) + II>. (84) 

3=1 



Zero-order approximation assumes that the fast dynamics in the closed system 
strongly couples the variables c, so that flows cannot influence this coupling. 

First-order approximationtak.es into account the shift of the invariant manifold 
by Sc. Equations for Newton's iterations have the same form (57) but instead 
of the vector field J they take into account the presence of the flow: 

[1 - P c ] (n + L' c 5c) = 0, P c Sc = 0, (85) 



where projector Pq corresponds to the unperturbed manifold. 

The first-order approximation means that fluxes change the coupling between 
the variables (concentrations). It is assumed that these new coupling is also 
set instantaneously (neglect of inertia). 

Remark. Various realizations of the first-order approximation in physical and 
chemical dynamics implement the viewpoint of an infinitely small chemical re- 
actor driven by the flow. In other words, this approximation is applicable in the 
Lagrangian system of coordinates (|Karlin, Gorban, Dukek fc Monnenmacher 



1998j ; Zmievskii, Kalin fc Dcvillc, 200C ). Transition to Eulerian coordinates is 



possible but the relations between concentrations and the flow will change its 
form. In a contrast, the more simplistic zero-order approximation is equally 
applicable in both the coordinate system, if it is valid. 
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11 Conclusion 



In this paper, we have presented the method for constructing the invariant 
manifolds for reducing systems of chemical kinetics. Our approach to com- 
putations of invariant manifolds of dissipative systems is close in spritit to 
the Kolmogorov-Arnold-Moser theory of invariant tori of Hamiltonian sys- 
tems ( [Arnold, 1963| , |1983j ): We also base our consideration on the Newton 
method instead of Taylor series expansions ( [Beyn fc Kless, 1998| ), and sys- 
tematically use duality structures. Recently, a version of an approach based 
on the invarinace equations was rediscovered by [Kazantzis (2000| ). He was 
solving the invariance equation by a Taylor series expansion. A counterpart of 
Taylor series expansions for constructing the slow invariant manifolds in the 
classical kinetic theory is the famous Chapman-Enskog method. The question 
of how this compares to iteration methods was studied extensively for certain 
classes of Grad moment equations QGorban fc Karlin, 1996a| ; [Karlin, Dukek & 



Nonncnmachcr, 1997a ; Karlin, 2000 ). 



The thermodynamic parameterization and the selfadjoint linearization arise 
in a natural way in the problem of finding slowest invariant manifolds for 
closed systems. This also leads to various applications in different approaches 
to reducing the description, in particular, to a thermodynamically consistent 
version of the intrinsic low-dimensional manifold, and to model kinetic equa- 
tions for lifting the reduced dynamics. Use of the thermodynamic projector 
makes it unnecessary global parameterizations of manifolds, and thus leads to 
computationally promising grid-based realizations. 

Invariant manifolds are constructed for closed space-independent chemical sys- 
tems. We also describe how to use these manifolds for modeling open and 
distributed systems. 
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Fig. 1. Images of the initial quasi-equilibrium manifold (bold line) and the first 
two corrections (solid normal lines) in the phase plane [01,03] for two-step catalytic 
reaction (63). Dashed lines are individual trajectories. 
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